Reduced chemodiversity suppresses rhizosphere microbiome functioning in the mono-cropped agroecosystems

Background Rhizodeposits regulate rhizosphere interactions, processes, nutrient and energy flow, and plant-microbe communication and thus play a vital role in maintaining soil and plant health. However, it remains unclear whether and how alteration in belowground carbon allocation and chemodiversity of rhizodeposits influences microbiome functioning in the rhizosphere ecosystems. To address this research gap, we investigated the relationship of rhizosphere carbon allocation and chemodiversity with microbiome biodiversity and functioning during peanut (Arachis hypogaea) continuous mono-cropping. After continuously labeling plants with 13CO2, we studied the chemodiversity and composition of rhizodeposits, along with the composition and diversity of active rhizosphere microbiome using metabolomic, amplicon, and shotgun metagenomic sequencing approaches based on DNA stable-isotope probing (DNA-SIP). Results Our results indicated that enrichment and depletion of rhizodeposits and active microbial taxa varied across plant growth stages and cropping durations. Specifically, a gradual decrease in the rhizosphere carbon allocation, chemodiversity, biodiversity and abundance of plant-beneficial taxa (such as Gemmatimonas, Streptomyces, Ramlibacter, and Lysobacter), and functional gene pathways (such as quorum sensing and biosynthesis of antibiotics) was observed with years of mono-cropping. We detected significant and strong correlations between rhizodeposits and rhizosphere microbiome biodiversity and functioning, though these were regulated by different ecological processes. For instance, rhizodeposits and active bacterial communities were mainly governed by deterministic and stochastic processes, respectively. Overall, the reduction in carbon deposition and chemodiversity during peanut continuous mono-cropping tended to suppress microbial biodiversity and its functions in the rhizosphere ecosystem. Conclusions Our results, for the first time, provide the evidence underlying the mechanism of rhizosphere microbiome malfunctioning in mono-cropped systems. Our study opens new avenues to deeply disentangle the complex plant-microbe interactions from the perspective of rhizodeposits chemodiversity and composition and will serve to guide future microbiome research for improving the functioning and services of soil ecosystems. Video abstract Supplementary Information The online version contains supplementary material available at 10.1186/s40168-022-01287-y.


Text S1. Calculation of the C allocation
The amount of photosynthesized 13 C accumulated in the plant organs (shoots and roots) and soils were calculated using Equation 1 [1]: 13 Camount(a)=[(atom 13 C%) l− (atom 13 C%) nl]a×TCa/100 (1) where "l" and "nl" are labeled and non-labeled samples, respectively; "a" refers to a plant organ (shoots and roots) or soil fraction; and "TCa" is the total carbon content (mg kg −1 ) of a.
Net assimilated 13 C was calculated using Equation 2: Net13Cassimilation(mg 13 C) = 13 Cshoots+ 13 Croot+ 13 Csoil (2) The photosynthesized 13 C distribution in plant and soil was calculated as follows: 4 domain signal acquisition. The signal to noise ratio and dynamic range were enhanced through accumulating 128 times domain FTICR transients.

Text S3. Isopycnic density gradient centrifugation and quantitative PCR
Total soil DNA was extracted from 0.5 g soil using a FastDNA SPIN Kit for Soil (approximately 340 μl each) were obtained for each sample, and the buoyant density of each fraction was determined by the refractive index using an AR200 digital hand-held refractometer (Reichert, Buffalo, NY, USA). The fractionated DNA was then precipitated using PEG 6000 for 2 h followed by centrifugation at 13,000 x g for 30 min. The DNA pellet was washed twice with 70% ethanol and finally dissolved in 30 μl TE buffer.
The 16S rRNA gene abundances from the DNA of each fraction isolated by the gradient centrifugation were quantified by quantitative PCR on a CFX96 Optical Real-Time Detection System (Bio-Rad, Hercules, CA, USA). The reaction was performed in a 20-μl mixture containing 10 μl SYBR Premix Ex Taq (Takara Biotech, Dalian, Liaoning, China), 0.5 μM of each primer, and 1 μl of DNA template. The primer pair 515F/907R was used to quantify the 16S rRNA gene, and the quantitative PCR conditions were as follows: 95 °C for 3 min, followed by 39 cycles of 95 °C for 30 s, 55 °C for 30 s, and 72 °C for 30 s [2]. Melting curve analysis was performed for every PCR analysis to confirm the specificity of the amplification products by continuously measuring the fluorescence as the temperature increased from 65 to 95 °C.
The relative abundance of 16S rRNA gene in the rhizosphere soils of unlabeled peanut plants ( 12 C-DNA) peaked at the buoyant density of 1.70 g mL −1 (light fraction) at both growth stages (W6 & W10); while that of labeled peanut plants ( 13 C-DNA) peaked at the buoyant density of 1.71 g mL −1 and 1.72 g mL −1 (heavy fraction) at W6 and W10, respectively (Fig. S1). The PCoA combined with the PERMANOVA based on the Bray-Curtis dissimilarity of bacterial communities indicated significant differences between the DNA groups (Fig. S1, Table S2).
These results clearly indicated that the 16S rRNA genes were successfully labeled, enriched, and separated by the isopycnic ultracentrifugation. The uncentrifuged DNA (raw DNA) had highest species richness, followed by 13 C light fractions, while 13 C light fractions, 12 C heavy fractions and 12 C light fractions showed no significant differences in the species richness (Fig. S1).
Raw sequence data were demultiplexed and quality filtered using the q2-demux plugin followed by denoising with DADA2 (via q2-dada2) [4], and the sequences that were not present in at least 2 samples were filtered out. After quality filtering and the removal of chimaeras, sequences were clustered into 2991 ASVs after rarefying to 15,333 sequences per sample (based on the sample with the minimum numbers of reads) [5]. The taxonomic assignment of representative sequences was performed using RDP classifier (http://rdp.cme.msu.edu/classifier/) with 80% confidence threshold [6]. A phylogenetic tree was first constructed using the Fasttree [7] after aligning the representative OTU sequences using PyNAST [8].

Text S5. Metagenome sequencing and data processing
The density-gradient centrifugation was conducted on DNA sample from each rhizosphere soil sample, and 13 C-DNA from the "heavy" gradient fractions was pooled to obtain sufficient 13 C-DNA for the metagenomics sequencing. Pairedend reads (2 × 150 bp) were generated on the Illumina Hiseq 4000 platform (Illumina Inc., San Diego, CA, USA). A total of 70.0-86.7 million raw reads per sample were yielded, with each sampling containing about 11.83 Gb raw data.
The resulting sequence reads were analyzed for quality control using a range of software programs: SeqPrep (https://github.com/jstjohn/SeqPrep) for tripping nonbiological bases in reads, such as primers or barcodes, and Sickle (https://github.com/najoshi/sickle) for filtering reads whose length after tripping was less than 50 bp and whose average quality score was less than 20.  [9], and those for which more than 90% of their length could be aligned to another gene with more than 95% identity (no gaps allowed) were removed to generate the non-redundant gene catalog.
The high-quality reads from each sample were aligned against the gene catalog by SOAPaligner (http://soap.genomics.org.cn/) using the criterion "identity > 95%." In each sample, the mapped reads of each gene were counted as the number of gene-mapped reads, and the gene relative abundance was calculated following a previous study. We aligned putative amino acid sequences, which were translated from the gene catalog, against the proteins/domains in KEGG databases (Release 79.0) using BLASTP (BLAST v.

= ∑ = ∑
Where ai is the relative abundance of gene i in sample S; Li represents the length of gene i; xi is the times at which gene i could be detected in sample S (the number of mapped reads); and bi is the copy number of gene i in the sequenced data from sample S.

Text S6. Detailed information about data analyses
The R package 'car' was used to analyze the normality of data distribution (Shapiro-Wilk test) and homogeneity of the variance (Levene's test). Then, the significant differences between treatments were determined using one-way analysis of variance (ANOVA), followed by Duncan's test for multiple comparisons (P < 0.05). However, if data were not normally distributed and/or variances were heterogeneous, the two-tailed Wilcoxon rank-sum test was performed to determine the statistical differences between treatments.
Three relational dendrograms based on molecular characteristics (MCD), potential biochemical transformations (TD), and transformation-weighted characteristics (TWCD) of the rhizodeposits were constructed [11]. Briefly, we of ultrahigh mass resolution variations between the identified rhizodeposits. We constructed a transformation network between rhizodeposits using the pairwise mass differences, while the correlations between rhizodeposits were calculated by picking the largest cluster of interconnected nodes (every node denotes an individual rhizodeposits molecule). Then, we measured the stepwise distance between each pair of the soil rhizodeposits. Afterwards, the TD was constructed using the UPGMA method, based on the standardized Euclidean distances. After combing the molecular characteristics and standardized transformation matrices, we constructed TWCD using the UPGMA method.